
# This code is hosted on http://code.google.com/p/lenthorp/
# Freely available for use in applications, but should NOT be modified
# Email all comments to lenthorpresearch@gmail.comimport RBFPricingEngine as rbfe

import RadialBasisFunctions as rbf
import numpy as np
import matplotlib.pyplot as plt
import BlackScholesImpliedVol as bs


from math import log
from math import exp
from math import sqrt
##simple Black-Scholes Euro option test

if __name__ == "__main__":

    spot   = 100
    strike = 100
    rate   = 0.05
    vol    = 0.35
    T      = 3.0
    n      = 51
    m      = 101
    sd     = 3

    xmin   = log(spot) -  sd * sqrt(T) * vol
    xmax   = log(spot) +  sd * sqrt(T) * vol
    dx     = (xmax - xmin)/(n - 1)
    ##build some nodes 1d
    nodes = []
    for i in range(n):
        x = xmin + dx * i
        nd= rbfe.Node(x)
        nodes.append(nd)
    ##set boundary nodes
    nodes[0].type  = 1
    nodes[-1].type = 2

    ##set boundary conditions
    boundCond    = dict()
    boundCond[1] = (lambda x,xc : rbf.MQ([x.centre],[xc.centre],sc,0), 0 )
    boundCond[2] = (lambda x,xc : rbf.MQ([x.centre],[xc.centre],sc,0), exp(xmax) - strike)


    ##time grid

    timeGrid = rbfe.TimeGrid(T,0.0,m)
    
    ##scaling coef
    sc = 4*dx
    
    ##create the equation
    ##2nd deriv term
    secDeriv       = lambda x,xc : rbf.MQ([x.centre],[xc.centre],sc,2)
    secDerivTerm   = (0.5*vol * vol , secDeriv )
    ##1st deriv term
    firstDeriv     = lambda x,xc : rbf.MQ([x.centre],[xc.centre],sc,1)
    firstDerivTerm =(rate - 0.5*vol*vol, firstDeriv)
    ##zeroth order term
    zeroDeriv      = lambda x,xc : rbf.MQ([x.centre],[xc.centre],sc,0)
    zeroDerivTerm  = (-rate,zeroDeriv)

    eqn = (secDerivTerm,firstDerivTerm,zeroDerivTerm)

    engine = rbfe.RBFPricingEngine(nodes,eqn,boundCond,timeGrid)
    v      = np.array([max(exp(x.centre) - strike,0) for x in nodes])
    engine.solveConstCoef(T,0,v,zeroDeriv)
    solution = engine.solution
    delta    = [ engine.calculateAtNode(firstDeriv,nd)/exp(nd.centre) for nd in nodes]
    gamma    = [ engine.calculateAtNode(secDeriv,nd)/(exp(2* nd.centre)) -
                 engine.calculateAtNode(firstDeriv,nd)/ (exp(2*nd.centre))   for nd in nodes]
    x        = [nd.centre for nd in nodes]

    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.set_ylabel('price')
    ax.set_xlabel('log S')
    ax.plot(x,solution)

    ax2 = ax.twinx()
    ax2.plot(x,delta)
    ax2.set_ylabel('delta')
    plt.show()

    spotNode = nodes[(n-1)/2]
    print "Solution at spot {0:0.3f}".format( engine.calculateAtNode(zeroDeriv,spotNode))
    fwd = exp(rate * T)*spot
    bsPrice = bs.BlackScholesPrice(fwd,strike,rate,T,vol) 
    print "Black Scholes price at spot {0:0.3f}".format( bsPrice )
    plt.show()
    

    
    
    

    

    
    
        
    
    